Example - Live-Cell Actin Filaments¶
This notebook is a demonstration of the super-resolution fluorescence imaging method, pysofi, on live-cell imaging data. The dataset is composed of 20 tiff videos in sequence ('Block1.tif' to 'Block20.tif'), and each video has 200 frames.
This notebook shows a standard workflow of pysofi to achieve high order Super-resolution Optical Fluctuation Imaging (SOFI) images. Briefly, we demonstrate a pipeline of
- Fourier interpolation
- Moments reconstrtuction or cumulants reconstruction
- Noise filtering 1
- Shrinking kernel deconvolution
- Noise filtering 2
- Local dynamic range compression (ldrc)
In the end, the processed image / video will be saved into a new tiff file with the colormap the user selects.
We also provide detailed introduction and demonstration of each step on separate notebooks, which can be found here.
import sys
sys.path.append('..')
import numpy as np
import matplotlib.pyplot as plt
from functions import visualization as v
from functions import (pysofi, filtering, masks, finterp)
import tifffile as tiff
%matplotlib inline
%load_ext autoreload
%autoreload 2
Load Image / Video (Tiff Files)¶
This sample dataset has 20 tiff files in an order from 1 to 20. First, we generate a list of filenames by their sequence, and load all videos into a dictionary dset, where keys are filenames in the list (e.g. 'Block1.tif') and values are PysofiData objects.
filenum = 20
filepath = '../SampleData'
filenames = ['Block' + str(i+1) + '.tif' for i in range(filenum)]
dset = {}
for filename in filenames:
dset[filename] = pysofi.PysofiData(filepath, filename)
If you don't get any message, the file is loaded successfully. We can check basic information of each data object including dimensions and frame numbers. For instance, for the third video ('Block3.tif'), we can get:
print("Information of '" + filenames[2] + "':\n")
print('Number of frames: ', dset[filenames[2]].n_frames)
print('Pixel numbers: ', dset[filenames[2]].xdim, '*', dset[filenames[2]].ydim)
We can also load one frame into an array with function get_frame(frame_num).
frame_im = dset[filenames[2]].get_frame(2)
plt.figure(figsize=(10,10))
plt.imshow(frame_im, cmap='pink')
1. Fourier Interpolation¶
This step performs Fourier interpolation to generate new pixels between physical pixels and reduce pixel sizes of each frame. The resulting TIFF file can be applied for SOFI processing to produce fSOFI results.
The user can carry out interpolation in a separate step and save the interpolated image / video for following pysofi steps, or skip this step and set finterp = True in moments or cumulants reconstructions.
We recommend to set the interpolation factor two times the reconstruction order number based on Nyquist sampling theorem. However, if your computation ability could not support too big interpolation factor, it is practical to use a smaller factor.
# Conduct Fourier interpolation on 'Block10' and save the interpolated image.
# You can skip this step if you do not need to explore or save interpolated images.
dset[filenames[9]].finterp_tiffstack(interp_num_lst=[4], frames=[0,10], save_option=True, return_option=False)
Instead of save the interpolated image stack as as '.tif' file, we can also return them in an array, and compare them with original images side-by-side using v.bokeh_visualization_mult(image_list, plot_title_list). As we can see below, the interpolated image have more pixels (1797 by 1797) comparing to the original image (450 by 450).
# Return the interpolated image (factor = 4) of the first frame of 'Block6.tif',
# and compare it with the original image.
# You can skip this step if you do not need to explore or return interpolated images.
finterp_im = dset[filenames[5]].finterp_tiffstack(interp_num_lst=[4], frames=[0,10], save_option=False, return_option=True)
original_im = dset[filenames[5]].get_frame(0)
# we can compare the images with and without fourier interpolation
plt.figure(figsize=(15,15))
plt.subplot(121)
plt.imshow(original_im,cmap='pink')
plt.title('original image')
plt.subplot(122)
plt.imshow(finterp_im[0][0],cmap='pink')
plt.title('Image with Fourier interpolation (interpolation factor of 4)')
As you can see, adding a Fourier interpolation step would increase the calculation time by a great amount. This applies to all following steps as well. Thus, based on the purpose and computation ability, we recommend the user to choose the appropriate interpolation factor. In each following step, we provide two options for the pysofi demonstration, and you can choose either one accordingly. Please choose the same option for all followinf steps.
- Option 1: pysofi on sixth-order moments-reconstruction series without Fourier interpolation.
- Option 2: pysofi on Fourier interpolated fourth-order moments-reconstruction series.
Option 1: sixth-order moments-reconstruction series without Fourier interpolation¶
2. Moments and cumulants reconstruction¶
Traditionally, the resolution enhancement of SOFI is manifested by the reduced width of the point spread function (PSF) in the cumulants-reconstructed SOFI image. However, one reason that limits the resolution enhancement of high-order SOFI is the positive and negative oscillations of the cumulants. The boundaries between the negative and positive regions generates “cusp artifacts” in high-order cumulant-reconstructed images. A practical (but mathematically nonrigorous) solution to avoiding cusp artifacts is to use moments reconstruction.
d.moment_image is a method of the PysofiData object that returns the moment-reconstructed images of a defined order. For more details on moment and cumulant reconstructions with pysofi, please see this notebook.
m_set = {}
for filename in filenames:
print('\n' + 'Calculating ' + filename + '...')
m_set[filename] = np.int64(dset[filename].moment_image(order=6, finterp=False))
#tiff.imwrite('moments_set6.tif', np.int64(np.around(m_set[filename])), dtype='int64',append=True)
# Plot M6 of 'Block6.tif'
v.bokeh_visualization(m_set[filenames[5]], save_option = False, imshow_same = True)
As we can see, high-order moments-reconstructions have a large dynamic range of pixel intensities. To overcome this problem, we conduct local dynamic range compression (ldrc) at step 6.
3. Noise Filtration 1¶
Moments-reconstruction is performed per block independently, each yield a single frame for further process. Due to the fact that the reconstructed image yielded from each individual block (1 to 20) is severely degraded by noise, we add a step of noise filtering along the time axis for each pixel independently by taking the convolution with a Gaussian kernel.
For more details on noise filtration using pysofi, please see this notebook.
# Option 1: Noise filtering on M6 seriers
nf = masks.gauss1d_mask(shape = (1,21), sigma = 2)
filtering.noise_filter1d(dset, m_set, nf, filtername = 'noise filter after moments-reconstruction')
# Show filtered array of block 6
print('\n' + filenames[5] + ' ' + dset[filenames[5]].filter_name)
print(dset[filenames[5]].filtered)
4. Shrinking Kernel Deconvolution¶
Deconvolution is then applied on each frame of the noise filtered movie to further enhance the spatial resolution. For the theory and demonstration of shrinking kernel deconvolution method, please see this notebook.
We can generate the original guess of the PSF (d.deconv_psf0 = masks.gauss2D_mask((401, 401), 2)), $\lambda$ (deconv_lambda = 1.5) and iteration number (deconv_iter = 20), and pass them to d.deconvsk.
# Option 1: Shrinking kernal deconvolution on filtered M6
deconv_psf0 = masks.gauss2d_mask((51, 51), 2)
deconv_psf0 = deconv_psf0 / np.max(deconv_psf0)
deconv_lambda = 1.5
deconv_iter = 20
deconv_set = {}
for filename in filenames:
print('\n' + filename)
im = dset[filename].filtered
deconv_set[filename] = dset[filename].deconvsk(deconv_psf0, im, deconv_lambda, deconv_iter)
#tiff.imwrite('deconvM6.tif', deconv_set[filename], append=True) # save tiff stack
5. Noise Filtration 2¶
Based on the same reasoning as before, a second noise filter (the same as noise filter 1) is carried out is applied per pixel along the time axis again to remove any possible irregular deconvolution artifacts. Note that the new filtered image is added to the PysofiData object as an attribute and will replace the previous filtered image.
# Option 1: Noise filtering on deconvoluted M6 image sereis
nf = masks.gauss1d_mask(shape = (1,21), sigma = 2)
filtering.noise_filter1d(dset,
deconv_set,
nf,
filtername = 'noise filter after deconvolution',
filenames=filenames)
6. Local dynamic range compression (ldrc)¶
High-order cumulants or moments reconstructions result-in images with a large dynamic range of pixel intensities. This ldrc algorithm compresses the dynamic range of these reconstructions with respect to a reference image (mask) while retaining resolution enhancement.The compression is performed locally in a small window that is scanned across the image. For details of the ldrc method, see Appendix 4 of this paper. We can pass the reference image mask_im = d.average_image(), scanning window size window_size = [25,25], image that needs compression (input_im = d.filtered) to the function. More information can be found here.
# Option 1: LDRC on the filtered deconvoluted M6 image sereis
window_size = [25, 25]
ldrc_set = {}
for filename in filenames:
print('\n' + 'Calculating ' + filename + '...')
mask_im = dset[filename].average_image()
ldrc_set[filename] = dset[filename].ldrc(window_size = window_size,
mask_im = mask_im,
input_im = dset[filename].filtered)
Visualization and Save Reuslts¶
After pysofi processing, we finally get super-resolution images stored in the dictionary ldrc_set. We can check the final images using Bokeh, save single frame into .png file or save the whole video into a .tif or .avi file. We can also enhance the contrast for better demonstration.
# Plot final image of block 4
v.bokeh_visualization(ldrc_set[filenames[5]])
# compare to average image:
v.bokeh_visualization_mult([dset[filenames[0]].average_image(), ldrc_set[filenames[0]]], ['Average', 'Pysofi'])
# Plot final image of block 10 with contrast enhancement
%autoreload 2
v.bokeh_visualization_rgba(v.enhance_contrast(ldrc_set[filenames[0]], display_contrast = 1.4))
# v.bokeh_visualization(ldrc_set[filenames[0]])
# save the final image stack as a .tif file
import tifffile as tiff
for filename in filenames:
im = ldrc_set[filename]
# stretch to cover the whole np.uint16 range
norm_im = (im - np.min(im)) / (np.max(im) - np.min(im)) * np.iinfo(np.uint16).max
tiff.imwrite(filepath + 'demo_pysofi.tif', np.uint16(norm_im), append=True)
# save the final image stack as a .avi file with enhance contrast
vid_array = [m_set[filename] for filename in filenames]
# For Windows, please use
# v.save_avi(vid_array, filename =filepath+'/demo_pysofi2.avi',os='windows')
v.save_avi(vid_array, filename = filepath + '/demo_pysofi.avi')
Option 2: Fourier interpolated fourth-order moments-reconstruction series.¶
2. Moments and cumulants reconstruction with finterp¶
Here, we can choose to conduct Fourier interpolation on each frame before the reconstruction by passing e.g. order = 4, finterp = True, interp_num = 8 to the function. Since Fourier interpolation would create almost 8*8 = 64 times more pixles in this case, this step might take some time. If your computation ability could not support too big interpolation factor, it is practical to use a smaller factor, for instance, interp_num = 2.
We recommend to check interpolated images with different factors by calling d.finterp_tiffstack with a list of interpolation numbers like interp_num_lst = [2,4,8] to check which interpolation order is desired, and pass only one factor for moments-reconstruction and further pysofi steps.
m_set_finterp = {}
for filename in filenames:
print('\n' + 'Calculating ' + filename + '...')
m_set_finterp[filename] = dset[filename].moment_image(order=4, finterp=True, interp_num=2)
#tiff.imwrite('finterp4.tif', np.int64(np.around(m_set_finterp[filename])), dtype='int64',append=True)
We recommend to use matplotlib.pyplot.imshow to show moments-reconstructions since Bokeh may crash Jupyter Notebook sometimes when the image size is too big. Only when we want to have an interactive image with zoom in/ out or hover tools, we would use Bokeh.
# Plot M4_finterp of 'Block11.tif' with palette 'Greens'
plt.figure(figsize=(8, 8))
plt.imshow(m_set_finterp[filenames[10]], cmap ='Greens')
3. Noise Filtration 1¶
# Option 2: Noise filtering on M4_finterp seriers
nf = masks.gauss1d_mask(shape = (1,21), sigma = 2)
filtering.noise_filter1d(dset, m_set_finterp, nf, filtername = 'noise filter after fourier interpolation and moments-reconstruction')
# Show filtered array of block 2
print('\n' + filenames[1] + ' ' + dset[filenames[1]].filter_name)
print(dset[filenames[1]].filtered)
4. Shrinking Kernel Deconvolution¶
Please change the size and sigma of the Gaussian kernel accordingly based on the Fourier interpolation factor. Note that when the image is too large (> 1500 * 1500 pixels), it is likely to come across an error when plotting it with Bokeh. In this case, we recommend to use matplotlib.pyplot.imshow instead.
# Option 2: Shrinking kernal deconvolution on filtered M4_finterp
deconv_psf0 = masks.gauss2d_mask((51*4, 51*4), 2*4)
deconv_psf0 = deconv_psf0 / np.max(deconv_psf0)
plt.imshow(deconv_psf0, cmap = 'pink')
deconv_lambda = 1.5
deconv_iter = 20
deconv_set ={}
for filename in filenames:
print('\n' + 'Calculating ' + filename + '...')
im = dset[filename].filtered
deconv_set[filename] = dset[filename].deconvsk(deconv_psf0, im, deconv_lambda, deconv_iter)
#tiff.imwrite('deconv.tif', deconv_set[filename], append=True) # save tiff stack
v.bokeh_visualization(deconv_set[filenames[1]][:,200:])
5. Noise Filtration¶
# Option 2: Noise filtering on deconvoluted image sereis
nf = masks.gauss1d_mask(shape = (1,21), sigma = 2)
filtering.noise_filter1d(dset,
deconv_set,
nf,
filtername = 'noise filter after deconvolution',
filenames=filenames)
v.bokeh_visualization(dset[filenames[3]].filtered[:,200:]) # filtered deconvoluted image of block 4
6. Local dynamic range compression (ldrc)¶
# Option 2: LDRC on the filtered deconvoluted M4-finterp image sereis
window_size = [25*4, 25*4]
ldrc_set_m4 = {}
for filename in filenames:
print('\n' + filename)
mask_im = dset[filename].average_image_with_finterp(4)
ldrc_set_m4[filename] = dset[filename].ldrc(window_size = window_size,
mask_im = mask_im,
input_im = dset[filename].filtered)
plt.figure(figsize=(10,10))
plt.imshow(ldrc_set_m4['Block1.tif'],cmap='pink')
